A homebuilt terahertz time-domain spectroscopy system was used as a THz radiation source. The 800-nm pulses with a repetition rate of 1 kHz from Ti:sapphire amplifier are focused into an LiNbO₃ crystal in tilted-pulse-front geometry to achieve velocity matching^[21,22] as shown in Fig. 1. THz pulses are generated through optical rectification, and the generated THz pulse has a time duration of 0.5 ps and peak electric field up to 150 kV/cm at the focus of two parabolic mirrors.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Schematic representation of experimental setup for exposure of cells to THz pulse radiation, displaying tilted-pulse-front THz pulse source. The THz beam is directed by flat metallic mirror to propagate to bottom surface of 12-well culture plate. M, PM, and represent flat mirror, paraboloidal mirror, and half-wave plate, respectively.

The THz wave intensity measured by electro–optic equivalent sampling in ZnTe crystal is shown in Fig. 2. The corresponding frequency domain spectrum obtained by applying fast Fourier transform is peaked at 0.7 THz with a full width at half maximum (FWHM) of 0.9 THz. The single pulse energy measured by opto-acoustic detector (Golay cell, Tydex) was close to . Thus, the average THz radiation power was evaluated to be approximately 1 mW, which yielded temperature increases less than 0.2 °C according to Pennes bio-heat transfer equation and the steady-state theoretical model.^[23] The broadband THz source allows the investigation of THz effects due to the high peak power up to 2 MW, when the low average power is maintained to eliminate thermal effects.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Waveform of (a) THz pulse intensity in time domain and (b) corresponding amplitude spectrum in frequency domain.

Three different types of eye cell lines were used (which were obtained from eye institute and school of optometry and ophthalmology, Tianjin medical university eye hospital), namely: human ARPE-19 retinal pigment epithelial cells (RPE for short), simian virus 40-transformed human corneal epithelial cells (HCE for short),^[24] and human MIO-M1 Müller cells (Müller for short). The cells were cultured in DMEM/F12 medium (Basalmedia, China) supplemented with 10% foetal bovine serum (FBS, Bovogen, Australia) and 1% penicillin/streptomycin solution (Gibco, USA) on culture dishes (NEST Biotechnology, China). Then, the cells were maintained in cell incubator (Thermo, USA) under standard conditions of 5% CO₂ and humidified atmosphere at 37 °C, while regularly changing the media and subculturing them onto fresh dishes when they are confluent by using trypsin (Gibco, USA).

The cells were seeded into 12-well cell culture plates 24 hours before THz exposure (NEST Biotechnology, China), with 1 ml of medium at a density of 1×10⁵ cells. For the control group, cells with the same passage were seeded adjacent to the exposed cells. To ensure that the cells were maintained in environments suitable for their long term survival, we used a custom-designed exposure chamber to keep the temperature at 36 °C–37 °C. The exposure chamber was made of insulating material, in which a circulator water bath device was used to control the temperature, assisted with an electronic thermometer (0.1 °C resolution) to measure the air temperature within the exposure chamber. A flat metallic mirror was used to direct THz radiation into the exposure chamber from below. The THz beam naturally diverged to a beam diameter that had approximately the size of the culture well and the corresponding average power density is nearly . Since the culture plate made of polystyrene has the low absorption of THz radiation,^[10] and the cells formed a monolayer adherent to the bottom, THz radiation could penetrate cells before being incident into the culture medium.

Cells with fresh culture medium were exposed to pulsed THz radiation for 6 hours, followed by incubation at 37 °C in 5% CO₂ atmosphere for 15 h without changing medium. Based on the calculation of gas flow rate, CO₂ was introduced into the exposure chamber with a duration of 12 s before THz irradiation, and the color of the medium did not change obviously due to the change of pH value. Owing to the large sample volumes required for RNA-seq, yields of four independent irradiation experimental replicates were pooled for post-processing.

The RNA-seq data are derived from sequencing of transcriptome and an alignment of the sequence to a reference genome, which yields highly reproducible results comparable to those of gene expression microarray.^[25] Total RNA was individually harvested from the irradiated cell cultures of three different cell types and their respective controls, using TRIzol reagent (Invitrogen, USA) according to the instructions of manufacturer. After assessing the purity, quantity and integrity of RNA samples, appropriate samples with an RNA integrity number (RIN) greater than 9.0 were subjected to library preparation. Messenger RNA (mRNA) was purified from total RNA with a value of as the input. Sequencing libraries were generated by using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, USA) following the recommendations. After cluster generation by using TruSeq PE Cluster Kit v3-cBot-HS (Illumia), the library preparations were sequenced on an Illumina Hiseq platform and 125 bp/150 bp paired-end reads were generated.

The raw data were purified to obtain clean data by removing low quality reads. The mapping of clean reads to the reference genome was performed by using STAR (v2.5.1b) based on the method of Maximal Mappable Prefix (MMP).^[26] The results of complete sequence alignment conditions as shown in Table 1, indicated that the uniquely mapping rates of all samples were greater than 93%. The fragments per kilobase million (FPKM) was subsequently calculated for estimating gene expression, based on the length of the gene and the reads count of clean data acquired by HTSeq v0.6.0 software.

Table 1.

Table 1.

Table 1. Complete sequence alignment conditions of samples. . Sample name Total reads Total mapping rate/% Uniquely mapping rate/% Multiple mapping rate/% RPE_Control 28115273 95.45 93.19 2.26 RPE_Exposed 25418095 96.33 94.02 2.31 HCE_Control 31357094 95.74 93.32 2.42 HCE_Exposed 29800525 95.74 93.32 2.44 Müller_Control 29658540 96.18 93.92 2.26 Müller_Exposed 28038821 96.02 93.69 2.33

Table 1. Complete sequence alignment conditions of samples. .

To investigate the long term influence of pulsed THz exposure on gene expression profile of cells, the analysis of differential gene expression is performed by using the edgeR R package (3.12.1) based on negative binomial distribution test.^[27] Benjamini and Hochberg false discovery rate (FDR) is used to eliminate the false positives, although it is unable to find all differentially expressed genes.^[15,28] Employing common significance criteria of FDR-adjusted statistical significance (p-values) , and expression fold change of two-fold down- or up-regulation at least, differentially expressed genes (DEGs) in three cell types are listed in Table 2.

Table 2.

Table 2.

Table 2. Differentially expressed genes. The list of genes identified by the RNA-seq when possessing differential expression with statistical significance (FDR-adjusted p-value ) of different cell types. . Cell type Ensembl ID Gene name Fold change (log2) FDR-adjusted p-value Gene type HCE ENSG00000138778 CENPE −2.163966507 0.018919368 protein-coding Müller nulla null null null null RPE ENSG00000165731 RET 9.267830464 2.32×10−8 protein-coding ENSG00000171094 ALK 4.05643875 1.05×10−6 protein-coding ENSG00000181143 MUC16 –5.671018261 5.24×10−6 protein-coding ENSG00000164692 4.673054049 5.62×10−6 protein-coding ENSG00000066468 FGFR2 5.012717306 7.27×10−6 protein-coding ENSG00000047936 ROS1 2.519821193 0.004482762 protein-coding a Null means none of differentially expressed genes selected in Müller cells.

Table 2. Differentially expressed genes. The list of genes identified by the RNA-seq when possessing differential expression with statistical significance (FDR-adjusted p-value ) of different cell types. .

For RPE cells, this analysis yields a total of 6 DEGs in the control versus exposed culture. In detail, the expression of RET, ALK, ROS1, FGFR2, and genes are up-regulated after THz irradiation, among which the largest fold change is recorded for RET. The RET, ALK, and ROS1 belong to the family of tyrosine kinase receptor, which play a role in cell growth, differentiation and migration.^[29,30] In addition, FGFR2 is a member of the fibroblast growth factor receptor family.^[31] Taken together, the observed up-regulations of these four genes may be suggestive of a higher level of sensitivity to the activation of a series of cellular processes. Indeed, it was previously shown that receptor tyrosine kinases are related to the development and progression of many types of cancers, while further investigations are needed to explore whether THz radiation can induce the RPE cells to transform into cancer cells.

Furthermore, there is only one statistically significant gene differentially expressed in HCE cells, named CENPE, which is down-regulated responded to the THz exposure. The CENPE is a kinesin-like motor protein, which is required for stable spindle microtubule capture at kinetochores during prometaphase of the cell cycle. The down-regulation of CENPE gene expression may be related to aneugenic effects observed in human fibroblasts exposed to THz radiation in some researches.^[10,32] While, in contrast to HCE cells, no gene in this significant criterion is affected by the THz exposure in Müller cells. In brief, the observations suggest that the influence of THz pulse on specific gene expression could last over 15 hours.

We examine the genes with less statistically significant evidence for their differential expression (p-value , not adjusted for multiple testing corrections, expression fold change ), which may provide more information about biological effects of THz radiation on a larger scale. To clearly visualize the global distribution of the DEGs out of over 23000 genes, volcano plots are presented in Fig. 3.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. Volcano plots of gene expression profiles in (a) RPE, (b) HCE, and (c) Müller cells. Fold change of expression magnitude is plotted versus level of statistical significance (transformed into a logarithmic scale). Note that red, green and gray color represent genes whose expression is considered up-regulated, down-regulated and unchanged, respectively. Besides, triangles represent the DEGs with FDR-adjusted p-value . Complete lists of DEGs, and corresponding fold change of expression and p-value are given in Data File 1.

Figure 3 shows a total of 297 DEGs detected with the cut-off criteria, 141 of which were up-regulated, 156 down-regulated in the control versus exposed cultures of RPE cells. For HCE cells, exposure to THz pulses changes the expression of 242 genes. The irradiation also affects the expression of 227 genes as compared with the control group in Müller cells. Overall, we find that THz radiation changes the transcription of approximately 1% of genes analyzed, which suggests that THz pulse may not have a strong or global influence on gene expression and its influence is only limited to a few genes. Meanwhile, none of the heat shock protein encoding genes are significantly up-regulated in exposed cells, which is in accordance with the scenario of the estimated temperature increasing as a result of THz exposure, demonstrating a minimum effect on gene expression caused by the THz-induced heating.

According to these data, now we come to analyze the overlaps in the plots of genes differentially expressed in three cell types. It is found that these genes are strongly different from each other in each of the cell types as shown in a Venn diagram (Fig. 4). While there is a little overlap between genes in every two of three cell types, no overlap exists among the gene expression profiles of all three cell types. We identify 18 DEGs that are common for two of three cell types, of which five genes are up-regulated and three genes are down-regulated due to the THz exposure. While, the other 10 genes display an opposite regulation when irradiated in different cell types. These considerable differences in gene expression profiles suggest that the influence of THz pulse is dependent on the type of cell. In addition, it may provide a hint toward the direct targets for the THz effects, which could exhibit the specificities of different cell types.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Venn diagram indicating down-regulated (green) and up-regulated (red) genes with p-value (not adjusted for multiple testing corrections) in these three cell types. Black color represents differentially expressed genes displaying opposite regulation in different cell types.

For further investigation on the specific biological functions affected by the THz radiation, the clusterProfiler R package (3.4.0) is used for gene ontology analysis based on Fisher’s exact test. Gene ontology (GO) represents an internationally standardized system for gene functional classification, covering three domains, i.e., cellular component, molecular function, and biological process. The functional classification of the DEGs with p-value reveals clear differences among the three cell types. Data in Fig. 5 are presented as the four GO terms with the highest level of statistical significance for each of the cell types.

Fig. 5.

Figure Option ViewDownloadNew Window

Fig. 5. Gene ontology enrichment analysis of differentially expressed genes. The first 4 GO terms with the highest level of statistical significance for three cell types are presented. Line graphs indicate the number of differentially expressed genes in specific GO terms, while bar graphs represent the level of statistical significance for enrichment. Complete lists of GO terms significantly enriched with FDR-adjusted p-value in three cell types are given in Data File 2.

Using the significance cutoff at the threshold value (FDR-adjusted p-value ), 53 functional GO terms are significantly enriched for the DEGs of RPE cells, among which five functional classes for molecular function are mainly related to the receptor protein kinase activity and motor activity, five GO terms belonging to cellular component domain are mainly in connection with extracellular matrix, and others about biological process are mainly linked to biological system development, differentiation and metabolic process. As for HCE cells, the GO analysis reveals statistically significant enrichment of DEGs in 14 GO terms concerning molecular function, which are mainly relevant to transmembrane transporter activity and ion channel activity. Only two enriched functional GO classes associated with biological process are identified for Müller cells, which are involved in metabolic process. The complete lists of GO terms significantly enriched with FDR-adjusted p-value in three cell types are given in Data File 2.

Alternative splicing (AS) is an important mechanism for regulating the expression of genes and the variable of protein, which results in a single gene coding for multiple proteins. The process of AS is regulated by a system of trans-acting proteins that bind to cis-acting sites on the primary transcript.^[33] To determine the AS events that are affected in response to THz pulses, rMATS(3.2.1) software is used to analyze the AS events based on likelihood-ratio test. The data shown in Fig. 6 illustrate the number of AS events influenced by THz radiation applying the common significant criteria (FDR adjusted p-value ). In RPE, HCE, and Müller cells, five, two, and six affected AS events are obtained, respectively, which are not uniformly distributed among the five types of AS events, but mainly concentrated in one or two types. The proportion of affected AS events is not more than 0.01%, which indicates that the effect of terahertz radiation on the AS event is specific. The complete lists of AS events significantly affected by FDR adjusted p-value in three cell types are given in Data File 3.

Fig. 6.

Figure Option ViewDownloadNew Window

Fig. 6. AS events analysis, showing the number of AS events influenced by THz radiation applying the common significant criteria (FDR adjusted p-value ) in three cell types. A3SS, A5SS, MXE, RI, and SE represent five types of AS events, alternative 3’ splice site (acceptor site), alternative 5’ splice site (donor site), mutually exclusive exons, retained intron and skipped exon, respectively. The complete lists of AS events affected by THz radiation with FDR adjusted p-value in three cell types are given in Data File 3.